Deep ecomorphological and genetic divergence in Steller's Jays (Cyanocitta stelleri, Aves: Corvidae)

Abstract The relationship between ecology and morphology is a cornerstone of evolutionary biology, and quantifying variation across environments can shed light on processes that give rise to biodiversity. Three morphotypes of the Steller's Jay (Cyanocitta stelleri) occupy different ecoregions in western North America, which vary in climate and landcover. These morphotypes (Coastal, Interior, Rocky Mountain) differ in size, plumage coloration, and head pattern. We sampled 1080 Steller's Jays from 68 populations (plus 11 outgroups) to address three main questions using data on morphology, plumage, genetics (mtDNA, microsatellites), and ecological niches: (1) How do phenotypic and genetic traits vary within and among populations, morphotypes, and ecoregions? (2) How do population‐level differences in Steller's Jays compare with other sister species pairs of North American birds? (3) What can we infer about the population history of Steller's Jays in relation to past climates, paleoecology, and niche evolution? We found substantial morphological, genetic, and ecological differentiation among morphotypes. The greatest genetic divergence separated Coastal and Interior morphotypes from the Rocky Mountain morphotype, which was associated with warmer, drier, and more open habitats. Microsatellites revealed additional structure between Coastal and Interior groups. The deep mtDNA split between Coastal/Interior and Rocky Mountain lineages of Steller's Jay (ND2 ~ 7.8%) is older than most North American avian sister species and dates to approximately 4.3 mya. Interior and Rocky Mountain morphotypes contact across a narrow zone with steep clines in traits and reduced gene flow. The distribution of the three morphotypes coincides with divergent varieties of ponderosa pine and Douglas fir. Species distribution models support multiple glacial refugia for Steller's Jays. Our integrative dataset combined with extensive geographic sampling provides compelling evidence for recognizing at least two species of Steller's Jay.

Morphological traits may exhibit both convergent and divergent patterns that are associated with ecology and behavior, and understanding these associations is important for studying drivers of diversification. For example, bird species living in coniferous forests show convergence in certain morphological traits (body mass, digital pads) that reflect adaptations to conifer needles (Korner-Nievergelt & Leisler, 2004), and similarities in avian morphological traits predict trophic level, dietary niche, and foraging behavior (Pigot et al., 2020). Conversely, local adaptation to different environments drives ecomorphological divergence and may promote reproductive isolation (Bertrand et al., 2016;Cicero & Koo, 2012;Ribeiro et al., 2014;Shakya et al., 2022). While ecological segregation is thought to drive rapid morphological change in some species (e.g., crossbills, Björklund et al., 2013), it is associated with morphological stasis in others (e.g., Empidonax flycatchers, Johnson & Cicero, 2002; Cinnyris Double-collared Sunbirds, McEntee et al., 2016McEntee et al., , 2021. Furthermore, different traits may vary in their evolutionary responses to environmental variation, as shown by distantly related clades of Myrmotherula antwrens that show similarities in body size but differences in body shape associated with ecology, behavior, and phylogenetic niche conservatism (Bravo et al., 2014). Studies that examine ecological and/or morphological divergence, especially in conjunction with molecular data, are therefore fundamental to elucidating patterns and processes of speciation. This is especially true for cryptic species, which may show subtle morphological differences among populations that occur in different environments (Bickford et al., 2007).
The western United States is well-suited to studies of ecomorphology because of its geological and paleoclimatic history and associated ecological diversity. Of the 72 ecoregions defined for the contiguous states, 36 (50%) occur from the Pacific Coast to the Rocky Mountains (west of longitude 105°; Olson et al., 2001). These various ecoregions collectively support many diverse species of both plants and animals. In a review of the biodiversity status of each state within the United States of America (Stein, 2002), eight of the 12 western continental states were ranked in the top 20 of all states in terms of species diversity, and 11 were ranked in the top 20 for endemism (Stein, 2002). This diversity is reflected in morphological, ecological, and genetic variation within and among closely related species, with new discoveries of locally adapted populations and divergent lineages uncovered regularly. For example, a new species of crossbill (Loxia sinesciurus) restricted to the South Hills and Albion Mountains of southern Idaho was described recently on the basis of morphological, genomic, behavioral, and ecological differences from other members of the Red Crossbill (L. curvirostra) complex (Benkman et al., 2009;Parchman et al., 2016). This species has evolved locally to specialize on the cone seeds of Rocky Mountain lodgepole pine (Pinus contorta latifolia) and presents a classic case of ecomorphological differentiation driving reproductive isolation from its close relatives. Similar integrative studies have identified ecomorphological divergence and genetic differentiation in other species of western birds, e.g., Dusky Grouse Dendragapus obscurus and Sooty Grouse D.
The Steller's Jay (Cyanocitta stelleri) is a common resident bird species of forests and woodlands in western North America and Central America. As many as 16 recognized subspecies are divided into three main groups (Clements et al., 2021;Walker et al., 2020): (1) the Coastal (stelleri) Group with four subspecies that range from southwestern Alaska through western British Columbia (including Haida Gwaii) and the United States west of the Cascade Mountains

T A X O N O M Y C L A S S I F I C A T I O N
Biodiversity ecology, Evolutionary ecology, Genetics, Phylogenetics, Population genetics, Taxonomy, Zoology and the Sierra Nevada; some authors (Browning, 2002) also recognize a fifth subspecies C. s. paralia in southwestern Oregon and northwestern California, but this taxon is synonymized with C. s. stelleri by others (Clements et al., 2021;Walker et al., 2020); (2) the Interior (diademata) Group with four subspecies that occur east of the Cascades and Sierra Nevada to the northern and southern Rocky Mountains, extending into northern Mexico; and (3) the Central American (coronata) Group with eight subspecies that range from the highlands of northern Mexico to northern Nicaragua. These three groups and associated subspecies are distinguished primarily by plumage coloration (breast and belly, back, head and crest), crest length, and facial patterning (extent and color of streaks on the forehead, which range from blue to white; presence and extent of a white patch or superciliary line above the eye). In addition, subspecies also differ in body size within and among groups, with a trend of larger birds in the north to smaller birds in the south (Walker et al., 2020). Steller's Jays occupy primarily coniferous and mixed coniferous-deciduous forests but are also found in pine-oak and oak woodland on exposures with cooler temperatures. Detailed phenotypic and ecological descriptions of Steller's Jay subspecies groups are provided by Walker et al. (2020).
Despite the extensive intraspecific variation, few studies have focused on ecomorphological and/or genetic variation in Steller's Jays.
The most complete morphological study to date was conducted by Brown (1963a), who examined 253 specimens (mostly adult males) from across the species' range to investigate ecogeographic factors influencing the crest as a specialized visual signal. Specific attention was paid to body size as measured by wing length, crest length and crest color, frontal (i.e., forehead) streaking, the contrast between the head and body coloration, plumage darkness and saturation, and length of the superciliary line. Results of this study showed strong geographic variation in these traits and associations with both climate and habitat. Brown (1963a) concluded that longer crests and more distinctive facial patterning (white frontal streaks and superciliary line), typical of the Interior Group, are associated with more open, drier forest habitats. Furthermore, he suggested that the frequency of conspecific encounters may be higher in these habitats, resulting in the selection of more conspicuous crests and facial patterns used as visual cues in aggressive encounters. Prior work by Brown (1960) showed that the crest is erected during aggression and serves as an important behavioral signal.
In another study, Bay (2002) compared morphological differentiation in relation to habitat in 156 Steller's Jays representing Rocky Mountain and California populations. In addition to finding significant between-group variability, results showed a strong association with habitat. For example, similar morphologies associated with feeding and locomotion were found in ecologically similar but geographically distant communities such as ponderosa pine (Pinus ponderosa) forest, suggesting an ecomorphological correspondence between sites.
Genetic data on Steller's Jays are currently limited to phylogenetic studies of the family Corvidae (Bonnacorso & Peterson, 2007;Erickson et al., 2005;Saunders & Edwards, 2000) and to a study of eight populations of Steller's Jay in the Pacific Northwest (Burg et al., 2005). The sister species of Steller's Jay is the Blue Jay (Cyanocitta cristata), and together they form a lineage that is part of a clade that includes the Pinyon Jay (Gymnorhinus cyanocephalus) and scrub-jays (Aphelocoma coerulescens, A.californica, A. woodhouseii) of North America; also related to this clade are Mexican to South American species of jays in the genera Calocitta, Cyanocorax, and Psilorhinus (Bonnacorso & Peterson, 2007;Erickson et al., 2005;Saunders & Edwards, 2000). Within the Steller's Jay, genetic data based on microsatellites for three northwestern subspecies (C. s. carolottae, C. s. stelleri, and S. s. annectens) found high levels of genetic variation among populations, with the highest divergence between the Haida Gwaii endemic C. s. carlottae and mainland subspecies C.
s. stelleri and C. s. annectens (Burg et al., 2005). Although C. s. stelleri and C. s. annectens were more similar genetically, these taxa also showed population structure with evidence for isolation by distance.
The high level of population differentiation found in northwestern Steller's Jays suggests rapid genetic differentiation following postglacial colonization in that region (Burg et al., 2005).
We sampled Steller's Jays extensively across their range in the western United States to investigate broader levels of genetic structure in association with a morphological and ecological variation.
Specifically, we focused on populations representing five subspecies that fall into three distinct morphotypes in the Coastal and Interior groups (Table 1). These populations represent the range of ecomorphological variation observed within our focal study area, which extends from the Pacific Coast to the Rocky Mountains. Because our specific goal was to examine phylogeography, population genetic structure, and phenotypic variation within and between these three morphotypes in relation to their environment, we did not sample across the entire range of subspecific and ecogeographic variation in Steller's Jays. Some subspecies and populations not included here have been the subject of prior study (e.g., C. s. carlottae and Pacific Northwest populations, Burg et al., 2005), while others are currently part of additional collaborative studies that extend south TA B L E 1 Primary trait characteristics distinguishing three morphotypes of Steller's Jays examined in this study.
into Mexico and Central America (Spellman et al., unpublished; McCormack et al., unpublished (Gavin, 2009;Shinneman et al., 2016) and experienced paleoclimatic events that have played an important role in divergence, speciation, and trait evolution in temperate birds (Friis et al., 2018;Johnson & Cicero, 2004;Lawson & Weir, 2014;Weir & Schluter, 2004). We examined variation in morphology, plumage, DNA, and environment to quantify patterns of differentiation within and among populations and morphotypes of Steller's Jays to address three main questions: (1) How do phenotypic and genetic traits vary within and among populations, morphotypes, and ecoregions? (2) How do population-level differences in Steller's Jays compare to other sister species pairs of North American birds? and (3) What can we infer about the population history of Steller's Jays in relation to past climates, paleoecology, and niche evolution? In addition, we evaluate these data in the context of Brown's (1963a) hypothesis that visual signaling traits (e.g., facial markings, crest length) have played an important role in Steller's Jay diversification.

| MATERIAL S AND ME THODS
We used an integrative approach to address the focal questions in our study. Data integral to understanding the results of our analyses are presented as figures or tables in the main text or in accompanying appendices (Figures A1-A14; Tables A1-A10). Additional supporting data are provided in supplementary materials (Tables S1-S7) and deposited in Dryad or Zenodo, along with the R code used to perform the analyses detailed below.

| Sampling
We sampled 1080 C. stelleri and assigned 1075 of those individuals to 68 populations from across the species' range in the western United States and Vancouver Island, British Columbia ( Figure 1); the remaining five individuals were from isolated locations with low sampling. Sampling of these 68 populations focused on five subspecies from the Coastal (C. s. stelleri, C. s. frontalis, C. s. carbonacea) and Interior (C. s. annectens, C. s. macrolopha) subspecies groups (Clements et al., 2021;Walker et al., 2020). In addition, we included 11 specimens of other C. stelleri subspecies or related species as outgroups: C. s. coronata from Mexico (1), C. s. ridgwayi from Mexico (1) and Guatemala (2), C. cristata (Blue Jay, 3), Aphelocoma californica (California Scrub-Jay, 1), Aphelocoma coerulescens (Florida Scrub-Jay, 1), and Gymnorhinus cyanocephalus (Pinyon Jay, 2). Of the 1091 total individuals sampled for this project (1080 ingroup samples of C. stelleri plus 11 outgroup samples), 1053 specimens are housed in the Museum of Vertebrate Zoology (MVZ), University of California, Berkeley. Data for another 33 specimens were obtained from specimens at the Burke Museum of Natural History and Culture, the University of Washington (n = 22, tissues), and the Museum of Southwestern Biology, University of New Mexico (n = 11, tissues plus vouchers). We also made use of five sequences sourced from GenBank. Details for the total dataset are given in Table S1.

| Morphological data collection and analyses
We used Fowler Ultra-Cal digital calipers to measure standard morphological characters on 1035 specimens (Table S2) prepared as study skins representing vouchers for the genetic samples; we limited our analyses to skins associated with tissue samples in order to associate the morphological data directly with the molecular data.
We only measured birds of known sex that had molted into at least their first-year (non-juvenile) plumage. Linear measurements (mm) included wing length, tail length, tarsus length, length of the middle toe, bill length, bill depth, and bill width (see Cicero, 1996:12 for details on methodology). We added the measurements for tarsus plus toe length to form a single character for analysis because the transition from tarsus to toe was not always clear on the specimens. We also measured crest length from the point of insertion of the longest crest feather to its tip. We were able to do this for all but a small (7%) number of birds, which had crests that were either too worn or molting heavily. Finally, we recorded the mass (grams) for every specimen from the data on the original specimen label.
We recorded the sex and age of each bird from the specimen label data. Sex was based on examination of reproductive organs during preservation and was scored as male, female, or unknown/ questionable. Age was recorded in two ways: (1) extent of skull ossification noted when the specimen was prepared; and (2) examination of tail shape, which can be used to distinguish first-year from older birds. Steller's Jays have a partial first prebasic (first year) molt, which involves the replacement of body but not wing and tail feathers (Walker et al., 2020). A small percentage (5%) of birds could not be aged easily by the tail shape, but they could be aged by ossification. Likewise, a few birds (less than 1%) could be aged by the tail but not by the skull. Age data were used to score each specimen as an adult (skull 100% ossified, adult tail) or immature (skull not fully ossified, first-year tail).
We recorded two plumage characters that are diagnostic of the different Steller's Jay morphotypes: (1) The color of frontal (forehead) streaking. Birds were scored as having either blue or white frontal streaks. Although Pacific Slope and Northwestern Interior birds have mostly blue streaking, close examination of the feathers revealed varying amounts of white mixed with blue. This contrasts with Rocky Mountain birds that have white streaking with no blue.
Thus, if the streaks contained any amount of blue, they were scored as being "blue." (2) The extent of the white superciliary line above the eye. We measured the length and width of the line using digital calipers and then multiplied width by length to estimate the eye line area. Birds with no white line were scored as "0" for both length and width. We recognize that a more quantitative approach to measuring plumage coloration (body plumage, eye-streak color) and patterning (size of frontal streaks) in Steller's Jays would be beneficial (Mason & Bowie, 2020), but these two characters provided a basic view of variation in traits that are important for distinguishing morphotypes.
To explore morphological variation among Steller's Jays by sex, age, morphotype, and population, we conducted a Principal Components Analysis (PCA) using the function princomp() from base R (R Core Team, 2020) and then implemented linear mixed models (LMMs) using the lme() function from the R package nlme v. 3.1.148 (Pinheiro et al., 2020). We included sex (male, female), age class (first year, adult), and morphotype (Coastal, Interior, Rocky Mountain) as fixed effects in all the models, and population as a random grouping effect because multiple individuals were typically sampled from the same population or locality. Using the R package MuMIn v 1.43.17 (Bartón, 2019), we calculated the marginal and conditional coefficients of determination (R 2 m and R 2 c , respectively) for each linear mixed model to assess whether including sampling locality as a random effect impacted the fit of the model (Nakagawa & Schielzeth, 2013).
For the LMMs, we omitted six populations from a putative contact F I G U R E 1 Sampling locations for Steller's Jays included in this study. Purple shading represents their range based on BirdLife International, edited to reflect known occurrences. Numbers correspond to the population numbers listed in Table S1. Circle colors correspond to the three distinct morphotypes (Table 1) present in the study area (yellow-Coastal; blue-Interior; purple-Rocky Mountain). Orange circles highlight a putative contact zone between Interior and Rocky Mountain populations based on the STRUCTURE analysis of microsatellites ( Figure 6). The circle sizes correspond roughly to the sample size at each location (see scale on the bottom left). The Steller's Jay plate on the bottom left (C. s. stelleri) represents the Coastal morphotype with a short crest, blue frontal streaks, and no supercilium. The Steller's Jay on the top right (C. s. annectens) represents the Interior morphotype with a short crest, blue frontal streaks, and white supercilium. The Steller's Jay on the bottom right (C. s. macrolopha) represents the Rocky Mountain morphotype with a long crest, white frontal streaks, and a white supercilium. Elevation is shown in grayscale with darker gray representing higher elevation. Illustrations are reproduced from https://birds ofthe world.org with permission from Lynx Edicions. zone ( Figure 1, Table S1) as they showed evidence of admixture from our ND2 haplotype network, STRUCTURE analysis, and phenotypic data (see below). We also generated box plots illustrating interquartile ranges for each morphological character and principal component axis for our fixed effect groupings (sex, age class, morphotype) and for each of the 68 populations. For our fixed effect groupings, we conducted a Tukey's honestly significant difference test with the R package agricolae v1. 3-5 (de Mendiburu, 2020) to compare the mean values of different groupings. Finally, we used the R package MASS v7.3-54 (Venables & Ripley, 2002) to conduct a Discriminant Function Analysis (DFA) with cross-validation using all continuous characters to quantify diagnosability among the three morphotypes (contact zone excluded).

| DNA laboratory methods
We extracted genomic DNA from frozen tissue samples using either a salt extraction protocol (Aljanabi & Martínez, 1997) or the Qiagen DNeasy Kit following the manufacturer's recommendations (Qiagen, Valencia, California). Following extraction, we used PCR to amplify the complete mitochondrial NADH Dehydrogenase subunit II (ND2) gene using primers L5204 and H6312 (Cicero & Johnson, 2001) with the following protocol: 94°C for 3 min, 35 cycles of 94°C denaturation for 30 s, 54°C annealing for 30 s, 72°C elongation for 45 s, and a final 10 min elongation at 72°C. We used 3 μl of each 10 μl PCR product for visualization on a 1% agarose gel stained with ethidium bromide. We then purified the remaining PCR products using 1 μl of 1:10 diluted ExoSAP solution incubated at 37°C for 30 min, followed by a denaturation step at 80°C for 15 min. We cycle-sequenced the purified products using Big Dye terminator chemistry (Applied Biosystems, Foster City, California), purified the cycle-sequencing products using Sephadex columns, and ran those on an ABI 3730 DNA Analyzer (Applied Biosystems). We assembled, checked, edited, and aligned all sequences using CodonCode Aligner version 4.2.7 (CodonCode Corporation).
We used 12 polymorphic microsatellite markers developed specifically for C. stelleri  to genotype 1075 of the 1080 ingroup samples that we also sequenced for ND2 (Table S4).
Of the five excluded individuals, four were analyzed for ND2 only (from GenBank sequences) and one (MVZ:Bird:179153) was dropped due to suspected contamination. All 12 loci were tetranucleotide repeats. We amplified the marker regions by PCR using the following reagents: 1× PCR Buffer (20 mM Tris HCl, 50 mM KCl, pH 8.4), 1.5 mM MgCl 2 , 10 mM of each dNTP, 10 mg/ml of BSA, 0.6 U of Taq polymerase, and approximately 10-30 ng genomic DNA. The thermocycling profile of the PCR reactions was 94°C for 2 min followed by 30-35 cycles of 94°C denaturation for 15 s, a locus-specific annealing temperature (53-57°C) for 15 s, and 72°C elongation for 15 s. We isolated 2 μl of the 10 μl PCR products for visualization on 1% agarose gels stained with ethidium bromide. We then mixed the PCR products with a GS-500 LIZ size standard (Applied Biosystems) and formamide, denatured the mixture at 94°C for 2 min, and loaded the mixture on an ABI 3730 DNA Analyzer. For 5 of the 12 markers, we cleaned the PCR product with T4 DNA polymerase (0.21 U) prior to loading to remove peaks of 3′ A nucleotide additions. We scored the resulting alleles for all 12 loci and 1075 individuals using GeneMapper 4.0 (Applied Biosystems).

| Phylogenetic analysis
We generated phylogenies from the ND2 alignment using both Bayesian and maximum likelihood methods for a subset of the Steller's Jays sampled for this study plus the 11 additional outgroup samples. Because of the large number of Steller's Jays in the entire dataset, we randomly selected two individuals from 10 populations representing each of the three morphotypes (total of n = 60 individuals). We determined the best-fitting model of nucleotide substitution with the programs PartitionFinder v2.1.1 (Lanfear et al., 2016) and PhyML (Guindon et al., 2010), and used the Bayesian information criterion to select from six models with each of the three codon positions as either linked or unlinked: JC, JC + I, HKY, HKY + I + G, GTR, GTR + G. The best performing model was HKY + I + G with separate, unlinked models inferred for each codon position partition. Using BEAST v2.4.5 (Bouckaert et al., 2014;Drummond & Rambaut, 2007), we then inferred a Bayesian phylogeny with a strict clock model of 2.1% divergence per million years (Weir & Schluter, 2008) that was linked across codon positions (see xml input file on Dryad for complete input settings). We ran three separate Bayesian BEAST analyses for 1 × 10 7 generations, from which we discarded the first 10% of generations as burn-in. We then combined the post-burn-in output of each run and generated a maximum clade credibility tree.
We also used RAxML v7.4.2 (Stamatakis, 2006) to infer a maximumlikelihood phylogeny. For RAxML, we used a GTR + G model of substitution, performed a rapid bootstrap analysis, and searched for the best-scoring maximum-likelihood tree within one program run ('-f a' option). Finally, we used the R package ape v5.0 (Paradis & Schliep, 2019) to calculate the uncorrected genetic distance among clades.

| Population genetic analyses
We generated a minimum spanning network of the ND2 haplotypes with PopART version 1.7 (Bandelt et al., 1999;Leigh & Bryant, 2015; see Table S3 for the Nexus input file). To reflect the geographic distribution of samples, we labeled each haplotype by morphotype as listed in Table S1.
We used a Bayesian approach to determine the level of population structure in the microsatellite data using the program STRUCTURE v2.3.4 (Hubisz et al., 2009;Pritchard et al., 2000).
To estimate the number of genetically distinct clusters (K), we used an admixture model with correlated allele frequencies and performed 10 independent runs of 5 × 10 5 MCMC iterations after a burn-in period of 1 × 10 5 iterations, using sampling locations as priors (LOCPRIOR). We assessed the likely number of K based on the inspection of plots ( Figure A1) for mean log likelihood of K (Ln P(K); Pritchard et al., 2000) and Delta K (ΔK; Evanno et al., 2005), both of which were obtained through Structure Harvester (Earl & vonHoldt, 2012). For the relevant runs, we used CLUMPP v1.1.2 (Jakobsson & Rosenberg, 2007) to generate mean membership coefficients for each individual sample using the "FullSearch" method (option 1).
To further visualize population clustering, we used the R package ADEGENET v2.1.1 (Jombart, 2008) to perform a discriminant analysis of principal components (DAPC, Jombart et al., 2010) on the microsatellite individual-genotype matrix. DAPC is a computationally fast, multivariate method for identifying the number of genetic clusters within a large dataset by maximizing between-group variation while minimizing within-group variation (Jombart et al., 2010).
To minimize overfitting, we performed an initial DAPC and evaluated the optimal number of principal components that would maximize the a-score; we used this number to select the number of principal components to retain in a subsequent re-analysis of the same dataset (Jombart & Collins, 2015). We then constructed a scatterplot of individuals by morphotype along the first two discriminant function axes to visualize separation based on microsatellites. We excluded populations from the putative contact zone, resulting in n = 1013.
We also used ARLEQUIN to estimate observed and expected heterozygosity (H O and H E respectively) and tested for departure from Hardy-Weinberg equilibrium at each of the 12 microsatellite loci (Excoffier & Lischer, 2010). We grouped samples by morphotype and excluded the contact zone populations. Results from the Hardy-Weinberg equilibrium test were corrected for multiple comparisons using the Holm-Bonferroni sequential correction (Holm, 1979).
We tested all populations and loci for linkage disequilibrium using GENEPOP v4.2 (Raymond & Rousset, 1995;Rousset, 2008). To account for the effects of population sampling on the number of alleles detected, we calculated the allelic richness (RS) for each locus in our populations using FSTAT version 2.9.3.2 (Goudet, 2002). We also calculated the number of private microsatellite alleles across all populations using the "Private Allele List" option in GenAlEx 6.5 (Peakall & Smouse, 2006. Finally, we used the 12 microsatellite loci to visualize genetic connectivity and variation among our sampled populations with the program EEMS (Estimated Effective Migration Surfaces, Petkova et al., 2016). We assigned individuals to one of 500 demes (i.e., equally spaced vertices) within the geographic extent of our study, and subsequently ran three separate Bayesian chains for 2 × 10 6 generations, of which we discarded the first 25% as burn-in. We examined stationarity and congruence among each of our runs and then visualized the output with the package rEEMSplots v0.9 (Petkova et al., 2016).

| Cline analysis
We used the R package HZAR v0.2.5 (Derryberry et al., 2014) to fit geographic clines for mtDNA, microsatellite, and morphological data across 11 populations from northern Idaho to southern Utah (41, 43, 49, 50-55, 60-61; Figure 1, Table A1) that showed divergence and evidence of admixture across the putative contact zone. This allowed us to estimate the center and width of clines using an MCMC algorithm to visualize changes in molecular and morphological variation and to determine potential mismatches between datasets. We calculated linear distances (rounded to the nearest 0.1 km) between the 11 populations using decimal coordinates and compressed all distance measurements into a single line on HZAR.
We fitted clines to PC1 scores of morphology, frequency of blue streaking, mtDNA haplotype frequencies, and Q scores obtained from STRUCTURE analyses of the microsatellite genotypes (see McEntee et al., 2016 and citations therein for use of Q scores in cline analyses). We fitted a total of 5 possible cline architectures that included a sigmoid curve with: (1) no exponential tails; (2) a right-flanking exponential tail; (3) a left-flanking exponential tail; (4) mirrored exponential tails; and (5) independent exponential tails on both sides. For each tail shape, we allowed the PMIN and PMAX (or "scaling") to be excluded (PMIN = 0 and PMAX = 1), fixed to the observed minimum and maximum, or free to vary. Thus, we tested a total of 15 possible models.
For each model, we ran a chain of 1 × 10 6 generations and discarded the first 5 × 10 5 generations as burn-in to optimize the covariance matrix used for an additional three independent MCMC chains (Derryberry et al., 2014). We ran the three subsequent chains for each model for a total of 9 × 10 6 generations and assessed convergence by inspecting the MCMC traces (Derryberry et al., 2014).
We selected the best model using the corrected Akaike Information Criterion (AIC c ).

| Comparison of ND2 divergence to sister species pairs
To place our study in a broader context, we extracted the raw percent divergence in ND2 sequences for a set of 30 North American avian species pairs (Table A2) and compared those values with the ND2 divergence between Coastal/Interior and Rocky Mountain Steller's Jays. We combined Coastal and Interior morphotypes for this analysis because together they form a separate mtDNA group from the Rocky Mountain populations (see below). We used the R package rentrez v1.2.3 (Winter, 2017) to query the National Center for Biotechnology Information (NCBI) nucleotide database and downloaded all available ND2 sequences for each species pair. After aligning each set of ND2 sequences using MUSCLE (Edgar, 2004) as called from the R package ape v5.4.1 (Paradis & Schliep, 2019), we calculated the minimum, mean, and maximum raw pairwise sequence divergence for each species pairwise comparison. We excluded sequences from individuals known to be in a contact zone between species (e.g., Aphelocoma californica and A. wollweberi, Gowen et al., 2014).

| Occurrence records
To examine how Steller's Jays are distributed across environmental space, we downloaded 1,032,960 specimen-based and observational occurrence records from 118 datasets accessed through the Global Biodiversity Information Facility (GBIF.org, 26 April 2019).
Using the R package sf (Pebesma, 2018), we subsequently excluded points that fell outside of the species' known geographic range, buffered by one degree, as well as those with a coordinate uncertainty (if provided) greater than 5 km. This resulted in 1,024,155 occurrence records.

| Sampling bias
Spatial sampling bias, which is pervasive in museum collections and observational databases (Beck et al., 2014), may occur when the relative sampling of environmental space is biased by the data collection protocol and does not represent the true environmental preferences of a species. Such bias is known to reduce the accuracy of species distribution models, especially when there is a lack of true absences (Boria et al., 2014;Fourcade et al., 2014;Kramer-Schadt et al., 2013;Phillips et al., 2009;Syfert et al., 2013). We took two approaches to lessen the impact of spatial sampling bias. First, we applied a proximity filter using the geoThin function in the R package enmSdm v0.5.3 (Smith, 2020) such that occurrences were no closer than 10 km from each other. When possible, we prioritized specimen-based occurrences over observational records. These filters resulted in 8262 occurrence records (Table S5). Second, we evaluated several anthropogenic predictors that could bias sampling, including distances to roads, populated places, urban areas, and protected areas. If such factors contribute to sampling bias in occurrence records, we wanted to account for the same bias in pseudo-absences by introducing it to the sampling of environmental space (Merow et al., 2016). We downloaded these datasets as vector data from www.natur alear thdata.com and found that species occurrences had the lowest mean distance to roads. Thus, we selected this anthropogenic predictor to represent sampling bias. After converting this distance to roads vector dataset to a grid format, we coded all grid cell values on a scale of 0 (minimal distance to roads) to 1 (maximal distance to roads). We then sampled 100,000 points in proportion to these grid cell values, such that coordinates near roads were more likely to be sampled. This set of points was used as pseudo-absences in subsequent species distribution modeling.

| Environmental predictors
We acquired climate data in the form of monthly precipitation and temperature grids from the CHELSA dataset v1.2 (Karger et al., 2017(Karger et al., , 2018. In addition to current climate, data were downloaded for the last glacial maximum (LGM, 21,000 years ago) under the PMIP3 NCAR-CCSM4 global circulation model, and for the period 2040-2060 under the CMIP5 NCAR-CCSM4 global circulation model (RCP4.5 scenario). All climate data were downloaded at a resolution of 30 arc-seconds and resampled to 60 arc-seconds (~ 2 km per grid cell).
For each time period, we generated 19 bioclimatic variables from the monthly CHELSA climate grids with the biovars function in the R package dismo v1.1-4 (Hijmans et al., 2017) and generated a set of 14 additional bioclimatic variables with the R package envirem v2.1 (Title & Bemmels, 2018). We used the R package raster v3.4-5 for general climate data manipulation (Hijmans, 2020). We also incorporated consensus landcover variables for the current climate period (landcover is not available for the LGM or future) in the form of nine landcover types (Tuanmu & Jetz, 2014); we excluded three landcover types (managed vegetation, urban, and open water) that were less relevant for our study. The full list of variables has been compiled in Table A3.

| Species distribution modeling
We first used the vifcor() function in the R package usdm v1.1-18 (Naimi et al., 2014) to reduce the collinearity in our set of predictor variables by assessing pairs of predictors that had a correlation coefficient greater than 0.8, and dropped the predictor with a greater variance inflation factor in each case. We then fit species distribution models with the maxnet algorithm , which is an implementation of the popular MaxEnt program as an infinitelyweighted logistic regression.
Overly parameter-rich models run the risk of being overfit to the data and of being less transferable to other time periods (Merow et al., 2014;Radosavljevic & Anderson, 2014;Warren & Seifert, 2011;Wright et al., 2014). Therefore, we further reduced the number of predictor variables and selected the optimal set of feature classes and smoothing parameters via a backward selection procedure that utilizes the corrected Akaike Information Criteria (AICc) as the evaluation criterion (Merow et al., 2013;Warren et al., 2014;Warren & Seifert, 2011). We started with all predictor variables (post initial reduction), identified the feature classes and smoothing parameters that led to the best AICc, and dropped the predictor variable that contributed the least to the model according to MaxEnt permutation importance. We repeated this until all included predictors had a permutation importance of at least 1%. For feature classes, we considered all combinations of linear, quadratic, and hinge features. We excluded product and threshold features because those were not found to be particularly important in the more recent versions of MaxEnt . We tested smoothing parameters from 1 to 12 with increasing increments of 0.5. This led to 161 parameter combinations. We ran these procedures for occurrences associated with each morphotype using a combination of existing and custom functions in R, with the help of packages enmSdm (Smith, 2020), dismo (Hijmans et al., 2017), maxnet , and raster (Hijmans, 2020).

| Discriminant function analysis
We performed a discriminant function analysis (DFA) with the environmental dataset using the MASS package v7.3-51 in R (Venables & Ripley, 2002) to determine the ability of climate and landcover variables to classify Steller's Jay occurrences into their respective morphotypes. We used the same reduced set of climate predictors selected for the SDM analyses as input variables and combined all forest-related classes (landcover classes 1-4) by summing the percent landcover to represent a single variable for closed habitat. We further used the DFA predictor loadings to identify the most influential variables in separating Steller's Jay morphotypes.

| Morphological variation
The PCA of morphological variation (Table A4) showed that the first axis loaded positively with body size (wing length, tail length, mass) and accounted for 43.42% of total variation. The second PCA axis, which accounted for 23.98% of the total variation, loaded positively with wing length, tail length, and crest length, and negatively with mass. Bill dimensions (length, width, or depth) did not load highly with either of the first two principal component axes.
Results from the linear mixed model analysis of morphological variation using principal component vectors (Table 2, Figure 2) showed statistically significant differences (p < .001) between sex and age groups, with males being larger than females ( Figure 2a) Table 2).
More details emerged when examining differences between each of the 68 populations analyzed in the linear mixed models (Figure 2d).
The most striking pattern was in the Coastal morphotype, which showed two parallel north-south clines in size (larger birds in the north) that corresponded to populations in the Coast Ranges (1-12) and Cascades-Sierra Nevada (13-38), respectively. Steller's Jays from populations representing the Interior morphotype (39-49) were relatively large and similar in size to birds from the more northern Coastal populations. Rocky Mountain birds (56-68) also showed a gradual trend toward smaller size in the south, but overall, they were similar morphologically to Steller's Jays from the central and southern part of the Coastal morphotype range. Morphological differences in individual traits and in PC2 all showed significant differences between at least two morphotypes (Table A5 and   Note: Population was included as a random grouping effect because multiple individuals were sampled from the same population or locality. p-Values below the .05 threshold are considered statistically significant. The effects of each factor are compared with a base model that describes variation among female, immature birds from coastal morphotype populations. The marginal (R 2 m ) and conditional (R 2 c ) coefficients of the linear mixed model are also shown.

TA B L E 2
Results from the linear mixed model of variation in the first principal component axis of morphological variation in Steller's Jays, which loads positively with overall body size.
Interior jays, and 80.0% of Rocky Mountain jays were classified correctly (Table 3). The 80%-96% success in classification indicates that morphological variation is partitioned geographically in Steller's Jays.

| Mitochondrial DNA variation
The ND2 phylogeny and haplotype networks showed congruent patterns of mitochondrial DNA variation within and between the three morphotypes ( Figure 5). In both the BEAST and RAxML analyses,  Population F I G U R E 3 Geographic variation in the frequency of (a) blue or white frontal streaks and (b) white supercilium among populations of Steller's Jays. Size of each circle is proportional to the number of individuals sampled from each population (as in Figure 1)

| Microsatellite variation and genetic connectivity
We successfully scored all 12 nuclear microsatellite loci for 1069 of the 1075 Steller's Jay individuals analyzed. One individual was successful for only 10 loci, and five individuals were successful for 11 loci (Table S1). All 12 loci were polymorphic across each of the 68 populations included in the analyses (  (Table S7). Therefore, we retained all loci in the dataset for further analyses.
Clustering analysis of the microsatellite data for all 68 populations using STRUCTURE showed an optimal value of K = 2 using the using Delta K (ΔK) method (Evanno et al., 2005; Table A7, Figure A1).

This is manifested as a sharp split of Coastal plus Interior versus
Rocky Mountain populations (Figure 6a), which is congruent with analyses of the ND2 sequence data. Membership coefficients from this analysis reveal a sharp transition from one group to the other across populations 50-55 in the putative contact zone.
While ΔK indicated an optimal K = 2 for the entire dataset, likelihood Ln P(K) values showed support for a higher K = 4 across all populations (Table A7, Figure A1) that separated the three morphotypes ( Figure 6b) and revealed some structure within Coastal populations (see STRUCTURE runs up to K = 7, Figure A12). All plots show more (c) DFA based on the reduced set of 13 climate predictors employed for the SDM analyses in addition to a landcover predictor representing closed habitat types. Each point represents an individual whose morphotype is defined according to subspecies and/or sampling locality. Contour lines encompass 95% of the data for each morphotype. The percent variance explained by each discriminant axis is provided in the axis labels.  F I G U R E 6 STRUCTURE results based on 12 nuclear microsatellite loci for 68 populations and 1075 individuals of Steller's Jay. Numbers correspond to populations in Figure 1 and Table S1, which were ordered to correspond with geography and ecomorphological phenotypes and assigned to clusters based on the Q scores from the STRUCTURE output. Bay (Grinnell, 1900;Grinnell & Miller, 1944). Within the Rocky Mountain morphotype, our dataset showed an optimal value of K = 2 using both methods. Although STRUCTURE plots ( Figure 6;   (Table 3).
The map of nuclear genetic diversity for Steller's Jays (Figure 7a) showed the highest levels in the central to southern Rocky Mountains The nuclear genetic diversity results are consistent with the private allele analysis (

F I G U R E 7
Estimated effective migration surfaces (EEMS) for Steller's Jays in the western United States, showing (a) nuclear genetic diversity and (b) effective migration rates; blue colors indicate higher genetic diversity and migration rates while orange colors indicate lower genetic diversity and migration rates, respectively. The polygon outlining our sampling localities represents the geographic extent of the analysis, within which populations were assigned to one of 500 demes represented by nodes. The size of each circle is proportional to the number of individuals assigned to a given node.

| Cline analyses
HZAR fits to the data showed strong and geographically coincident clines ( Figure 8) in ND2 haplotype frequencies, blue frontal streak frequencies, and microsatellite Q scores between populations 50-55 that span 421 km. The microsatellite data showed the steepest cline, which was offset slightly to the south compared with the other two traits; this break occurred between populations 51 and 53, which are located 129.7 km apart. Although morphological PC1 scores also showed a general cline in this region, the curve was much shallower and occurred over a larger geographic distance to the south (populations 53-60, 215 km). A comparison of the different HZAR models and associated values is given in Table A8.

| Environmental variation
Collinearity reduction of climatic variables decreased the set of predictors from 33 to 12. When landcover variables were included, the initial set of 42 predictors (33 climates plus 9 landcovers) was reduced to 22. A list of predictor variables is given in Table A9. Feature classes and smoothing parameters for MaxEnt, optimized via AICc model selection, can be found in Table A10.
Species distribution models for the three morphotypes ( Figure 9) showed strong geographic patterning. The model based on the current climate and landcover closely matched present-day distributions and separated the morphotypes into corresponding Coastal,  Figure 1 and are plotted above the x-axis to show the location and distance of symbols on each cline. Note that the cline for morphology PC1 is on a different scale and axis than the other three traits.    (Ashton, 2002), and may result from both phenotypic plasticity and genetically-based adaptations (Stillwell, 2010). Bay (2002) showed that morphological differentiation in Steller's Jays is correlated with habitat and appears to be adaptive to local environments. Studies of other phenotypically variable taxa (e.g., Junco hyemalis) have shown that local adaptation driven by environment can promote rapid differentiation (Friis et al., 2018). Our finding of differentiation in central Coastal California provides evidence of the historical isolation of the Steller's Jay subspecies carbonacea described over a century ago (Fisher, 1902;Grinnell, 1900). This subspecies was named on the basis of size and color characters intermediate between stelleri and frontalis, and its distribution was described originally as coastal Oregon and California from the Columbia River south to Monterey County (west of the Cascades; Grinnell, 1900;Fisher, 1902). Later, its range was restricted to central Coastal California based on examination of newly collected specimens (including those in fresh plumage) from northwestern California (Maillard, 1922). Maillard (1922) noted that the central coastal area occupied by carbonacea "was supposed at one time to be either an island or a group of islands not widely separated," and hypothesized that Steller's Jays spread from the interior mountains of California toward the coast where they diverged in isolation as an insular form. Support for this hypothesis comes from other taxa that show similar genetic and/or phenotypic divergence in central coastal California. For example, the songbird subspecies Junco hyemalis pinosus that breeds in this area is differentiated across its genome from other junco taxa (Friis et al., 2022), and this geographic region is a hotspot of mammalian (Davis et al., 2008) and herpetological (Rissler et al., 2006) diversification as well as plant endemism (Stebbins & Major, 1965). Extensive sampling and whole genome sequencing across diverse taxa (including Steller's Jays) for the California Conservation Genomics Project (Shaffer et al., 2022) will shed further light on lineage diversification in central coastal California as well as other regions and habitats in the state.

F I G U R E 9
The divergence between Coastal, Interior, and Rocky Mountain populations of Steller's Jays in mtDNA and microsatellites is geographically coincident with splits between other montane avian taxa. Barrowclough et al. (2004)

| Ecomorphological associations
Results from this study integrating morphology and environment on a broad geographic scale, combined with evidence that birds in similar habitats (e.g., ponderosa pine) have similar morphologies in different geographic regions (Bay, 2002), provide strong evidence of ecomorphological associations in Steller's Jays. Discriminant function analysis and species distribution models for Steller's Jay occurrence records showed that the three morphotypes occur in distinct ecoregions based on both climate and landcover, and that 92%-97% of the records were correctly classified to morphotypes according to these ecological variables. Furthermore, occurrence density plots based on climate and percent landcover showed the Rocky Mountain morphotype to be more common in areas of lower moisture (Climatic Moisture Index), higher minimum temperature of the warmest month, higher monthly and annual temperature ranges varying light environments appear to be a major factor driving avian plumage evolution, with selection favoring either conspicuousness or crypsis (Goméz & Théry, 2007;McNaught & Owens, 2002;Shultz & Burns, 2013). In Steller's Jays, it is possible that brighter and more conspicuous facial features used for visual signaling in open habitats (Rocky Mountain morphotype) may be counterbalanced by selection for darker features (e.g., blue frontal streaks, no superciliary line) that are more cryptic in lower light, closed coastal forest habitats (Coastal morphotype). Plumage traits of individual Coastal Steller's Jays have been found to influence the rate of extra-pair parentage and the proximity of territories to forest edge (Overeem et al., 2014).
While studies have shown that selection in different habitats can be an important driver of phenotypic divergence, the outcome may be clade, environment, or sex-specific (Cicero et al., 2020;Cornuault et al., 2015;Mason & Bowie, 2020;Medina et al., 2017).
Population density and social behavior also may influence visual signaling in Steller's Jays. In addition to habitat openness, Brown (1963a) noted that conspecific encounters should vary as a function of population density and intraspecific competition.
Data from the North American Breeding Bird Survey Sauer et al., 2017) show that the highest density of Steller's Jays occurs in the range of the Coastal morphotype, with notably lower densities in the distribution of both Interior and Rocky Mountain morphotypes.
Relative abundance during the breeding season varies widely across the species' range, from a high of 46.4 birds per survey route in the Sierra Nevada to a low of 0.1 bird per route in Wyoming (Walker et al., 2020). Personal field experience (CC) corroborates this regional difference in population density-during the collection of specimens for this study, Coastal Steller's Jays were encountered routinely in higher numbers and in social groups compared with further east where the birds were less common and often solitary or in pairs. Studies in coastal California have shown that Steller's Jays exhibit dominance hierarchies that reflect social rank maintained by frequent interactions (Brown, 1963b), and that socially monogamous pairs have extensively overlapping home ranges with neighboring birds congregating at food resources (Kalinowski et al., 2015).
Furthermore, Coastal Steller's Jays cache food items in a social context, whereby birds travel larger distances to cache when other birds are present (Kalinowski et al., 2015). Although comparable data are lacking for Interior and Rocky Mountain populations, one might expect Steller's Jays to encounter other individuals less frequently in areas of low population density, especially if they do not occur in social groups. The interaction between population density, habitat openness, behavioral encounters, and visual signaling traits across morphotypes of Steller's Jays requires further study.
Finally, the role of vocalizations for signaling is a critical component of understanding lineage divergence and ecomorphological associations. Song complexity is significantly associated with habitat openness in New World sparrows, which also show stronger phylogenetic signal in behavioral traits (song structure, vocal duets) compared with morphological traits (Cicero et al., 2020). By contrast, visual and vocal signals reportedly evolve independently, with no habitat association, across the diverse family of tanager species (Mason et al., 2014). Steller's Jays are notoriously loud and have a wide variety of vocalizations and mimetic calls that are used in both territorial defense and predator response (Billings et al., 2017;Brown, 1964;Hope, 1980;Walker et al., 2020). Vocal variation among Steller's Jay populations or morphotypes has not been stud-

| Population history
This study presents mtDNA and microsatellite data that shed light on the population history of Steller's Jays in the western United States. Divergence estimates based on mtDNA place the split of the Steller's Jay and Blue Jay in the late Miocene or early Pliocene, approximately 6.36 million years ago (95% highest posterior density = 5.20-7.46 mya). Steller's and Blue jays occupy different bioclimatic regimes in western and eastern North America, respectively, with Blue Jays breeding primarily in more humid hardwood habitats east of the Rocky Mountains including in the Great Plains (with some recent westward expansion, Smith et al., 2020). The estimated timing of divergence of these two sister species coincides with a trend toward progressive drying and cooling in the Great Plains that continued to the Pleistocene (Burke et al., 2018;Frye & Leonard, 1957).
The Great Plains is well-known as an area of divergence and secondary contact between western and eastern lineages of birds and other taxa (Reding et al., 2021;Rising, 1983;Swenson & Howard, 2005).
Although these two species of jays are known to hybridize (Walker et al., 2020), such events are probably infrequent (Williams & Wheat, 1971).
Within the Steller's Jay, Coastal/Interior and Rocky Mountain lineages show an exceptionally deep split for avian species (7.8% in mtDNA) that dates to approximately 4.3 mya, which places their divergence in the mid-Pliocene during a period of warmer conditions (Burke et al., 2018). This split also is reflected in the microsatellite clusters. Interestingly, the distribution of these two lineages is coincident with divergent Coastal and Rocky Mountain varieties of ponderosa pine (Latta & Mitton, 1999;  Mountains (Gugger et al., 2010). Ponderosa pine, which is more arid-adapted, is hypothesized to have occurred in southern coastal, southern interior highland, and Sierra Madre refugia (Roberts & Hamaan, 2015). Microsatellite data from Steller's Jays in the Pacific Northwest (coastal Washington through British Columbia to southeastern Alaska) revealed rapid divergence and postglacial expansion from at least one coastal refugium (Burg et al., 2005). Our  Baird (1854:118) recognized Rocky Mountain populations of Steller's Jay as a distinct species in his original taxonomic description of Cyanocitta macrolopha. Approximately 20 years later, Coues (1871) discussed his experience with the "Long-crested Jay" macrolopha and commented that this species and the Steller's Jay are "so much alike that they might be considered as one species." That comment may have been the rationale for lumping the two taxa in the first edition of the checklist to North American Birds, which recognized frontalis and macrolopha as distinct subspecies (American Ornithologists' Union [AOU], 1886). Although that treatment has stood for over a century, we strongly support re-elevating Rocky Mountain populations to species status.

| Species limits and taxonomic recommendations
Species delimitation requires an integrative taxonomy approach with strong geographic sampling that includes putative contact zones . In this study, we combined multiple lines of evi- Concept . In addition to their deep genetic divergence, these populations are morphologically and ecologically distinct with evidence of limited gene flow between them. A comparison of mtDNA sequence divergence between pairs of North American avian species showed that Coastal/Interior and Rocky Mountain Steller's Jays rank among the most divergent, and cline analysis showed a sharp break from southwestern Wyoming and southeastern Idaho to northern Utah. The steepest cline was observed in the microsatellite data over a distance of approximately 130 km (between populations 51 and 53), which matches the average cline width for avian hybrid zones in North America (Slager et al., 2020). Finally, the niche modeling data showed that Coastal/Interior and Rocky Mountain morphotypes occupy different bioclimatic regimes, have a preference for different landcover densities (closed and open, respectively), and experienced separate histories in glacial refugia.
The splitting of species requires consideration of alternate English names, as outlined by the Guidelines for English Bird Names published by the American Ornithological Society's Committee on Classification and Nomenclature of North American Birds (https:// ameri canor nitho logy.org/nacc/guide lines -for-engli sh-bird-names).
Specifically, the guidelines state that for "true phylogenetic daughter species formerly treated as a single parental species, the usual policy is to create a new name for each daughter species. This practice is designed to prevent confusion in the literature as to what taxonomic entity the parental name…references." Prior to 1957, polytypic species of North American birds had English names assigned to subspecies rather than to the species itself (Strickland, 2017). In the case of the Steller's Jay, two alternate names are well-established in the literature: Blue-fronted Jay (AOU, 1886) or Blue-fronted Steller Jay (Grinnell & Miller, 1944) for the Coastal subspecies frontalis, and Long-crested Jay (Baird, 1854;AOU, 1886) or Long-crested Steller Jay (Grinnell & Miller, 1944) for the Rocky Mountain subspecies macrolopha. These names are descriptive of the blue frontal streaking that distinguishes the Coastal/Interior morphotype, and the long crest typical of the Rocky Mountain morphotype. Although Whitefronted Jay is another descriptive name for Rocky Mountain populations and presents a nice contrast with Blue-fronted Jay, that name has not been established in the literature. In keeping with the English Name Guidelines reference above, and for simplification of the names, we recommend the scientific and English names Bluefronted Jay (Cyanocitta stelleri) and Long-crested Jay (Cyanocitta macrolopha) for Coastal/Interior and Rocky Mountain morphotypes, respectively.
In addition to the populations studied here, our mtDNA data revealed that populations from southern Mexico and Guatemala were strongly divergent (uncorrected ND2 sequence divergence of 3.76%) from those to the north in the Rocky Mountains-a value higher than the median (2.15%) for other North American avian species that we compared. Steller's Jays from those southern populations are phenotypically distinct (Brown, 1963a), most notably with a blue rather than black crest. Because our current study did not focus on these southern populations and included only a few samples as outgroups, we cannot recommend further taxonomic action at this time. More detailed study with dense sampling is needed to assess genetic differences, timing of divergence, and potential zones of intergradation between Rocky Mountain, Mexican, and Central American populations. Although these populations are part of the Rocky Mountain clade, we suspect that such a study will provide evidence for additional splitting in future.

AUTH O R CO NTR I B UTI O N S
Carla Cicero: Conceptualization (lead); data curation (equal); funding acquisition (equal); investigation (equal); methodology (equal); project administration (lead); supervision (lead); visualization (equal); writing -original draft (lead); writing -review and editing (equal).  (Table S5), and the museum staff and collectors who helped track down specimens for GenBank accessions (Table A2) This paper is dedicated to Ned K. Johnson, whose knowledge and love of western North American birds stimulated this study.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflict of interest.

O PEN R E S E A RCH BA D G E S
This article has earned an Open Data badge for making publicly available the digitally-shareable data necessary to reproduce the reported results. The data is available at https://doi.org/10.6078/ D14Q5N.

DATA AVA I L A B I L I T Y S TAT E M E N T
DNA sequences are deposited in GenBank, accessions OM689560-OM690612, OM817568-OM817600 (Table S1)

S U PP O RTI N G I N FO R M ATI O N
Additional supporting information can be found online in the Supporting Information section at the end of this article.
How to cite this article: Cicero, C., Mason, N. A., Oong, Z Note: Population numbers correspond to those in Figure 1 and Note: In each model, population was included as a random grouping effect because multiple individuals were sampled from the same population or locality. p-Values that fall below the .05 threshold are considered statistically significant and shown in bold. The effects of each factor are compared with a base model that describes variation among female, juvenile birds from Coastal morphotype populations. The marginal (R 2 m ) and conditional (R 2 c ) coefficients of the models also are shown.   Population F I G U R E A 1 2 STRUCTURE plots of microsatellite data (K = 2 through K = 7) for all 68 populations of Steller's Jays sampled in this study. Numbers correspond to populations in Figure 1 and Table S1, which were ordered to correspond with geography and ecomorphological phenotypes and assigned to clusters based on the Q scores from the STRUCTURE output. Colors of the different genetic clusters are the same as in Figure 6.

F I G U R E A 1 3
Absolute values of DFA loadings for 13 climate and 1 landcover variables. See Table A8 for further description of the predictors.